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Abstract. In this paper, we develop a formalism in order to incorporate the contribution of internal gravity waves to the transport 
of angular momentum and chemicals over long time-scales in stars. We show that the development of a double peaked shear 
layer acts as a filter for waves, and how the asymmetry of this filter produces momentum extraction from the core when it is 
rotating faster than the surface. Using only this filtered flux, it is possible to follow the contribution of internal waves over long 
£vq ] (evolutionary) time-scales. 

(^-^ . We then present the evolution of the internal rotation profile using this formalism for stars which are spun down via magnetic 

torquing. We show that waves tend to slow down the core, creating a "slow" front that may then propagate from the core to the 

l/^ , surface. Further spin down of the surface leads to the formation of a new front. Finally we show how this momentum transport 

reduces rotational mixing in a 1.2 M e , Z - 0.02 model, leading to a surface lithium abundance in agreement with observations 

"*^h ' in the Hyades. 
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j3 I 1. Introduction one has to pay attention to their possible interactions and inter- 

• • i dependence 
. Stellar models are getting more and more sophisticated. In or- 
"O 1 der to explain detailed observed features of stars in various These mechanisms affect evolutionary tracks, lifetimes, 

K i places of the Hertzsprung-Russell diagram, modern stellar evo- surface abundances, chemical yields etc. Their impact on the 

13 . -, ,. . . • j j ■ . , , u • i stellar structure and evolution arises mainly from the redistribu- 

lution codes must indeed incorporate several complex physical 3 

, • , f . e j . u . i jjj rrii tion of chemical elements that they cause, and which is related 
processes which are often referred to as non-standard . The . .... 



main ones are : 



(except in the case of atomic diffusion) to the induced redistri- 
bution of angular momentum inside the stars. During the last 

- Atomic diffusion (gravitational settling, thermal diffusion, fifteen years, many studies have been devoted to describe the 
radiative forces); evolution of the angular momentum distribution because this 

- Large scale mixing due to rotation (turbulence, meridional pattern governs the extent of rotation-induced mixing in stellar 
circulation); interiors. For models in which the internal rotation law evolves 

- Convective overshooting; under the effects of meridional circulation, shear mixing, hor- 

- Internal gravity waves (IGW); izontal turbulence, mass loss, contraction and expansion (i.e., 

- Magnetic fields. neglecting the transport by IGW and by magnetic fields) the 

main results can be summarized as follows : 

All these processes are not necessarily present at the same 

time everywhere in the HR diagram. Additionally they act on „ , , , , , • , . , . , . , 

3 b 3 3 - On the one hand, models which take into account the hy- 

very different characteristic timescales in stars of various initial , ■ ■ , • j j u j-uj 

3 drodynamical processes induced by rotation as described 

masses and evolutionary stages. In order to provide a coherent fey Zahn (1992)> Maedej . (1995)> Talon & Zahn (199?) 

picture of stellar evolution, one must thus understand why some ■ , „ ^ , /ln no\ ft- i • 
F 3 and Maeder & Zahn (1998) are very successful in explain- 
dominate in certain stellar types and not in others. Furthermore, . , , „ ,. „ . . ,. , , . , 
3r ing the observed anomalies in stars which do not have 

Send offprint requests to: Suzanne Talon extended envelopes. Among the many successes of these 
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models one can quote the reproduction of the left side of 
the Li dip in field and open cluster stars, of the constancy of 
the CNO abundances within the Li dip and of the evolution 
of Li abundance in sub-giants (Charbonnel & Talon 1999, 
Palacios et al. 2003, Pasquini et al.2004). For more massive 
stars these models explain for example the observed He and 
N enrichment in main sequence O- and early B-type stars, 
in OB super-giants and in A-type SMC super-giants. They 
also account for the observed variations of the Wolf-Rayet 
star populations as well as the fractions of type Ib/Ic su- 
pernovae with respect to type II SN at various metallicities 
(see Maeder & Meynet 2000, Meynet & Maeder 2005 and 
references therein). 
- On the other hand, the same input physics fails to repro- 
duce some of the most constraining observed features in 
low-mass stars which have deep convection envelope. First 
of all, meridional circulation and shear turbulence alone 
are not efficient enough to enforce the flat solar rotation 
profile measured by helioseismology (Brown et al. 1989; 
Kosovichev et al. 1997; Couvidat et al. 2003; Pinsonneault 
et al. 1989; Chaboyer et al. 1995 1 ; Matias & Zahn 1998). 
Additionally the rise of the Li abundance on the red side of 
the dip cannot be explained if one assumes that momentum 
is transported only by the wind-driven meridional circula- 
tion in those main-sequence stars which are efficiently spun 
down via magnetic torquing (Talon & Charbonnel 1998). 
Last but not least, the current rotating models are insuffi- 
cient to explain the anomalies in C and N isotopes observed 
in field and cluster red giant stars (Palacios et al. 2005). 

The successes and difficulties described above have re- 
vealed the occurrence of an additional process that participates 
in the transport of angular momentum in relatively low-mass 
stars which have extended convective envelopes and are spun 
down via magnetic braking in their early evolution (Talon & 
Charbonnel 1998). This process must act in conjunction with 
meridional circulation, turbulence and atomic diffusion. As of 
now, two candidates have received some attention, namely in- 
ternal gravity waves and magnetic fields. 

Internal gravity waves have initially been invoked as a 
source of mixing for chemicals (Press 1981, Garcia Lopez & 
Spruit 1991, Schatzman 1993, Montalban 1994, Montalban & 
Schatzman 1996, 2000, Young et al. 2003). Ando (1986) stud- 
ied the transport of momentum associated with standing, grav- 
ity waves. He showed how momentum redistribution by these 
waves may increase the surface velocity to induce episodic 
mass-loss in Be stars. Goldreich & Nicholson (1989) used them 
later in order to explain the evolution of the velocity of binary 
stars, producing synchronization that proceeds from the surface 
to the core. Traveling internal gravity waves have since been 
invoked as an important process in the redistribution of angu- 
lar momentum in single stars spun down by magnetic torquing 
(Schatzman 1993, Kumar & Quataert 1997, Zahn, Talon, & 
Matias 1997). 

1 Note that the Yale group computed the evolution of angular mo- 
mentum using a simplified description of the action of meridional cir- 
culation which was considered as a diffusive rather than as an advec- 
tive process 



In a previous series of papers, we examined the genera- 
tion of internal gravity waves by the surface convection zone 
of stars with various masses and metallicities. We found that 
these waves, which are able to extract angular momentum from 
the deep solar interior (Talon, Kumar, & Zahn 2002, hereafter 
TKZ), have a very peculiar mass (or more precisely effec- 
tive temperature) dependence and could possibly dominate the 
transport of angular momentum in stars with deep enough con- 
vective envelope (Talon & Charbonnel 2003, 2004). We sug- 
gested that such a dependence could lead to a coherent picture 
of rotational mixing in stars of all masses at various evolution 
phases. It could for example explain simultaneously the cold 
side of the Li dip as well as the solar rotation profile and the 
existence of fast rotating horizontal branch stars (see Talon & 
Charbonnel 2004 for more details). 

An important characteristic of internal waves is that, un- 
less they are damped, they conserve their angular momentum 
even when their local frequency is modified by Doppler shift- 
ing. This should be kept in mind throughout the reading of this 
paper. For a comprehensive review of gravity wave properties, 
we suggest Bretherton (1969) and Zahn et al. (1997) for an ap- 
plication to the stellar, spherical case. 

The other mechanism that has been invoked to enforce the 
Sun's flat rotation profile requires a pre-existing fossil magnetic 
field (Charbonneau & Mac Gregor 1993; Barnes, Charbonneau, 
& Mac Gregor 1999). However no mass dependence is ex- 
pected in that case. This is in contradiction with the Li dip con- 
straint (although in this case this feature could be explained in 
purely solid body rotation by the combined action of solid body 
meridional circulation and radiative forces; see Charbonneau & 
Michaud 1988) and more importantly with the large rotational 
velocities measured in some horizontal branch stars. 

Let us also mention calculations of massive star models by 
Maeder & Meynet (2003) in which the Taylor-Spruit dynamo 
(Spruit 1999, 2002), which is thought to be the most efficient 
dynamo mechanism in the radiative region, is included. When 
considered as is, this mechanism induces almost solid body ro- 
tation, efficiently reducing the extent of rotational mixing. It 
also leads to stellar models whose properties are closer to those 
of standard models and thus, in lesser agreement with stellar 
observations. When reviewing the efficiency of this mechanism 
by taking into account energy conservation, it is found that 
the dynamo is not as efficient as first thought. In fact, calcula- 
tions for a 15 M Q star model show that adding the contribution 
of magnetic field only slightly modifies the results obtained 
with the purely hydrodynamic models (Maeder & Meynet 2004 
and Maeder, private communication). Let us recall that for 
such massive stars with no convective envelope internal gravity 
waves can not be produced efficiently and thus the present ro- 
tating models would remain unchanged in our global picture. 2 

All these results and constraints obtained for stars cover- 
ing a large fraction of the HR diagram strongly incite further 
studies of the combined effects of rotation and internal gravity 
waves in stars where these waves are known to be efficiently 

2 The potential effect of internal gravity waves produced by the con- 
vective stellar core should be modest, although it remains to be quan- 
tified. 
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Fig. 1. Schematic view of the physical processes and their interactions involved in momentum redistribution by IGW. Complete 
calculations take into account all processes enclosed in the dashed box. 



produced. Although they do not definitively disgrace magnetic 
fields which certainly play a role in the complete picture, they 
bring sufficient arguments to justify complete tests of a fully 
hydrodynamical model before including MHD effects. 

This is what we intend to do in the present series of pa- 
pers where we will investigate the combined and sometimes 
conflicting effects of rotation, internal waves and atomic dif- 
fusion in low- and intermediate-mass stars at various stages of 
their evolution. We will present the first fully hydrodynamical 
stellar models that include self-consistently the chemical and 
momentum redistribution by these three mechanisms. 

In the present paper we discuss the formalism that we will 
use in the hydrodynamical models which will be presented here 
and in forthcoming studies. A schematic view of all physical in- 
gredients included and their interactions is shown in Fig.^ In 
§ 2, we first describe the prescription we follow for wave gener- 
ation or excitation. Then in § 3, we come back to the properties 
of the shear layer oscillation (or SLO) which builds up below 
the surface convection zone of the star and on the momentum 
extraction by the waves in the radiative region. We explain how 
an excited wave spectrum combined with the action of the SLO 
may be replaced by a filtered spectrum and a diffusion coeffi- 
cient, as illustrated in Fig. ^ General equations for the trans- 
port of angular momentum and chemicals are written in the 
global scheme (§ 4 and 5), which takes into account momen- 
tum deposition beyond the shear layer and momentum trans- 
port by meridional circulation and shear turbulence. Then we 
recall some properties of the transport by waves using a static 
model (§ 6). Finally we present the first results for a 1.2 M , 
Z = 0.02 star computed within the complete physical frame- 
work (§ 7) before concluding (§ 8). 

2. Internal Gravity Wave Generation 

Gravity waves are produced, among other things, by the injec- 
tion of kinetic energy from a turbulent region to an adjacent 
stable region. This is observed for example at the border of 



clouds in the earth's atmosphere (Townsend 1965) and also in 
laboratory experiments (Townsend 1958). It should also occur 
in stars. This was already illustrated in the early (2D) numeri- 
cal simulations of convection including penetration by Hurlburt 
et al. (1986, see also Hurlburt et al. 1994, Andersen 1994, 
Nordlund et al. 1996, Kiraga et al. 2000, Rogers & Glatzmeier 
2005). There are two ways to excite those waves: 

- convective overshooting in a stable region; 

- bulk excitation, similar to that of the solar pressure waves. 

Excitation by overshooting is probably the most difficult 
to evaluate analytically. The first attempt was made by Garcia 
Lopez & Spruit (1991) by assuming that the pressure perturba- 
tion produced by turbulent eddies at the radiative/convective 
boundary is equal to the wave pressure perturbation. In this 
model, stochastic eddies of a given size are considered to con- 
tribute to the excitation of a whole wavelength spectrum. As 
formulated, this model rests on the assumption of homoge- 
neous turbulence, although it could be modified to incorporate 
the presence of downdrafts observed in numerical simulations 
of convection. Another estimate for this process has been made 
by Fritts, Vadas, & Andreassen (1998), in a study aimed at es- 
timating the residual circulation induced by latitude dependent 
dissipation in the tachocline. However, they concentrate on the 
small wavelength waves that dissipate close to the convection 
zone and their mechanism does not take into account the com- 
bination of small scale eddies in order to produce low degree 
waves. Such waves are essential in order to influence the inner 
layer on an evolutionary time-scale (TKZ). Let us also mention 
here that the results presented in this study depend on the latitu- 
dinal differential rotation and cannot (as of now) be generalized 
to stars other than the Sun. 

Kiraga et al. (2003) tried to evaluate the validity of the 
Garcia Lopez & Spruit approach by comparing the predictions 
of this model to a simulation of penetrative convection. The 
peak wave spectrum produced by the simulation was similar 
in amplitude to that of this parametric model. However, in the 
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Fig. 2. Penetration of IGW into the convection zone of a ZAMS 
1 M Q model for various degrees I. The arrows indicate the 
depth where the amplitude is reduced by a factor of 2. Also 
shown is the local convective flux given by the mixing length 
theory. 



flux just below the convection zone, v is the convective veloc- 
ity, L is the radial size of an energy bearing turbulent eddy, 
Ti w L/v is the characteristic convective time, and is the 
radial size of the largest eddy at depth r with characteristic fre- 
quency of co or greater (h m = Lminjl, (2wt7,)~ 3 / 2 }). The radial 
wave number is related to the horizontal wave number by 

where A^ 2 is the Brunt-Vaisala frequency. In the convection 
zone, the mode is evanescent and the penetration depth varies 
as VF(7TT) 3 . Figure |2 compares this distance for various 
modes with the convective energy that is locally available. In 
the outer part of the model, the luminosity is carried almost 
uniquely by convection. 

The momentum flux per unit frequency Tj is then related 
to the energy flux by 

Tj{m,e,o)) = -r E (t,cS) (3) 
to 

(Goldreich & Nicholson 1989, Zahn et al. 1997). 



simulation modes where excited over a much broader range of 
frequencies and wavelengths. It is not clear whether this is re- 
lated to the bi-dimensional nature of the simulation of to short- 
comings of the model. 

Internal gravity waves can also be excited in the convection 
zone itself. In that region, modes are evanescent and their am- 
plitude is proportional to exp y— Jdr^j, where k r is the radial 
wave number. Press (1981) used the formalism of Goldreich 
& Keeley (1977) to show that Reynolds-stress and buoyancy 
may excite gravity waves with a rather large amplitude at the 
bottom of the convection zone. Goldreich & Kumar (1990) 
and Goldreich, Murray, & Kumar (1994, GMK) completed the 
Goldreich & Keeley formalism in order to apply it to solar 
p-modes. Their model quite successfully reproduces the solar 
spectral energy input rate distribution, provided one free pa- 
rameter which describes the geometry of turbulent eddies is 
calibrated. In that case, driving is dominated by entropy fluctu- 
ations. Balmforth (1992) made a similar study, using a some- 
what different formalism. Subject to the calibration of a free 
parameter, he is also able to reproduce the spectral energy dis- 
tribution; however, it is the Reynolds-stress that is the main 
source of driving. 

In the present study, we follow Kumar & Quataert (1997) 
and apply the GMK formalism to traveling internal gravity 
waves. This was also used by Kumar, Talon, & Zahn (1999), 
TKZ, Talon & Charbonnel (2003, 2004). The energy flux per 
unit frequency Te is then 

™»>-£/*?[(i) 1+ « + »(*f 

x«p[-/4«< + l)/2^] TT ^ B5 , (1) 

where and [({( + l)] 1 ^ 2 ^;, are the radial and horizontal dis- 
placement wave-functions which are normalized to unit energy 



3. Evolution of Angular Momentum 

3.1. Wave Mean-Flow Interaction and the SLO 

It is now well established that the dissipation of internal gravity 
waves in a differentially rotating region leads to an increase in 
the local differential rotation. In stellar models, this leads to the 
formation of a narrow (~ 1-2% in radius) doubled peak oscil- 
lating shear layer adjacent to the convection zone where they 
are produced (Gough & Mclntyre 1998, Ringot 1998, Kumar, 
Talon, & Zahn 1999). We shall refer to this process as "shear 
layer oscillation" or SLO. 

In a simple two waves model, Kim & MacGregor (2001) 
examined the behavior of that layer, which depends on the ratio 
between viscosity and wave flux. If the viscosity is large, a sta- 
tionary solution may be found. There then exists a bifurcation 
to an oscillatory behavior of the shear as viscosity is reduced. 
Further reduction leads to the appearance of chaos. The deter- 
mination of the stellar viscosity is thus of key importance for 
the structure of this layer. 

If only radiative viscosity is considered, one observes the 
formation of a very steep and narrow shear layer, and chaotic 
behavior is expected. However, as the local shear increases, it 
may lead to the appearance of a shear instability, which will 
enhance the local viscosity compared to its microscopic value. 
This process actually self-regulates the wave-mean flow inter- 
action. Indeed, a larger wave flux leads to a larger differential 
rotation which in turns acts as to increase the local viscosity. 

Several calculations have been performed, using various 
prescriptions for the turbulent viscosity. The first case-study 
should consist of using simply the radiative viscosity. It has 
the major advantage that it can be derived from first principles 
only. However, this viscosity is so low compared to the wave 

3 This is the theoretical dependence. Damping can be enhanced by 
turbulence, which would reduce the surface amplitude of the mode 
(see e.g. Andersen, 1996). 
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Fig. 3. Behavior of the shear layer as a function of the diminishing viscosity for a 1.3 M Q model with a turbulent viscosity given 
by Eq. @. From left to right, N = 8, N = 7.9 and N = 7.85. {top) Time evolution of the rotation velocity at two different depths. 
(Full line) Q(n) (dashed line) £2(>2) with r\ = r cz - 0.003Hp and ri = r c: - O.OlSHp. (bottom) Phase space diagram. 



flux that it rapidly leads to a slow layer that is stopped and even 
begins to rotate backward. To prevent that, we will consider a 
"turbulent" viscosity proportional to the radiative viscosity 

V = lO^Y-ad. (4) 

A regular oscillating layer is formed for N = & (a similar form 
was used by TKZ). For a slightly larger (N = 9) viscosity, we 
no longer obtain an oscillation. This corresponds to the station- 
ary solution found by Kim & MacGregor (2003). Here, we do 
not find a stationary profile due to the use of different bound- 
ary conditions (Neumann rather than Dirichlet); our solution 
slowly decays to a state of solid body rotation. For smaller 
(N = 7.9 and N = 7.85) viscosities, there is a transition to 
chaos. These results generalize the Kim & MacGregor (2001) 
results to the case where a complete wave spectrum is consid- 
ered. Period doubling is not readily identified in this simula- 
tion. Figure|3]summarizes our findings. 

The turbulent viscosities used so far are not realistic. They 
vary only slowly with depth, and remain large even far from re- 
gions where there exists a physical mechanism to produce this 
turbulence. We must thus find a reliable prescription in order to 
estimate the magnitude of turbulence on physical grounds. 

The structure of the shear region points to the shear insta- 
bility as an important source of turbulence. Generally, its mag- 
nitude depends on the local shear rate and on the efficiency of 
buoyancy which acts as a restoring force. When taking into ac- 
count radiative losses (Townsend 1958, Dudis 1974, Lignieres 
et al. 1999) and the effect of horizontal turbulence on the sta- 
bilizing effect of mean molecular weight gradients (Talon & 
Zahn 1997), the viscous turbulence may be approximated by 

8 Ri ait (rdn/dr) 2 
Vv 5N 2 T /(K + D h ) + NpD h ' 

This formulation is based on the local shear rate. In a 
rapidly varying profile, such as the one we get in the shear 



layer, it leads to large gradients in the turbulent diffusion co- 
efficient, and is probably not realistic. Indeed, vertical turbu- 
lence can develop over vertical scales similar to the pressure 
scale height which is greater than the scale of the variations of 
the shear layer (less than 1% in radius, see e.g. Fig.0] while 
in low mass main sequence stars the pressure scale height just 
below the convection zone ranges from about ~ 10% in ra- 
dius in a 1 M star to ~ 3% in a 1.4 M star). In practice, we 
thus used a convolution of the local shear turbulence given by 
Eq. with a Gaussian of width 0.2 H P (H P is the pressure 
scale height) in order to obtain a more realistic turbulent coeffi- 
cient 4 . Turbulence is also averaged over a complete oscillation 
cycle. 

In this framework, the magnitude of the shear layer is self- 
regulated. An increase in the wave flux leads to an increased 
differential rotation and thus, an increase of the shear turbu- 
lence. 

Calculations are shown here for a 1.3 M , Z = 0.02 model. 
Figure |4] shows the evolution of the shear layer as well as 
the self-consistent turbulent viscosity induced in the model. 
In this model, thermal diffusion is rather large (see Talon & 
Charbonnel 2003). This explains the large value for the tur- 
bulent coefficient, and also its extreme thinness. In lower mass 
stars, while just below the convection zone the viscosity is large 
enough to lead to a regular oscillating layer, it then decreases 
rapidly in the interior, leading to a complex layer structure that 
has a chaotic behavior. However, this has no particular impact 
on long term momentum extraction. 

This turbulence is generated by wave momentum depo- 
sition due to radiative damping, and is related to turbulence 
produced by large amplitude waves where they are dissipated 
(Press 1981, Canuto 2002, Young et al. 2003). It is quite dif- 

4 The Gaussian we use here is quite thin. Far from a sharp shear 
layer, the convolution we are using would not alter the turbulent vis- 
cosity from that obtained with Eq. 
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Fig. 4. (left) Shear layer evolution (successive profiles are sep- 
arated by 1 year) and (right) average turbulent viscosity (as de- 
scribed in § 13.11 . The dashed line corresponds to an analytical 
fit (cf. Eq.EJ. 



ferent in essence from the turbulent diffusion generally asso- 
ciated with waves, either due to wave breaking (Press 1981, 
Garcia Lopez & Spruit 1991) or to irreversible second order 
motions (Press 1981, Schatzman 1993, Montalban 1994). This 
issue will be discussed further in §|5] 

3.2. Momentum Deposition Beyond the Shear Layer 

The superficial shear layer acts as a filter on internal gravity 
waves (Gough & Mclntyre 1998). Indeed, as internal waves 
travel across first a "rapid" and then a "slow" layer, prograde 
and retrograde waves respectively are preferentially damped. 
However, some of the power remains in the lowest order waves, 
and over long time-scales (of the order of evolutionary time- 
scales), momentum redistribution below the shear layer cannot 
be neglected (TKZ). Figure [5] illustrates wave characteristics 
below the shear layer for a given differential rotation, with the 
convection zone rotating slower than the radiative zone. The 
two top rows compare the amplitude of retrograde and prograde 
waves. As shown in the third panel, the retrograde waves are 
somewhat less damped then their prograde counterparts, and 
this leads to a net luminosity of momentum, illustrated in the 
bottom panel. As explained in TKZ this is because the under- 
lying differential rotation produces a prograde layer which, on 
average, is always larger than the retrograde layer. 

It is the low frequency, low degree waves that give the 
largest contribution to momentum redistribution in the inte- 
rior; high degree waves are damped closer to the convection 
zone (as damping oc g (( + l) 3 ^ 2 ) and thus have an intrinsically 
small amplitude below the shear layer while high frequency 
waves experience less differential damping (oc l/cr(r) 4 = 
1/ [to - m (Q(r) - £2 CZ )] 4 ), crucial to produce a net momentum 
deposition (see Fig.^Jl. See §|4]for details on damping. 

In order to built evolutionary stellar models, it is not possi- 
ble to follow in detail the behavior of the oscillating shear layer, 
which occurs on time-scales of years or tens of years. However, 




Fig. 5. Wave characteristics below the shear layer for a differ- 
ential rotation of 6Q. = 0.1 //Hz over 5% in radius in a 1.3 M 
star. First and second row. Wave amplitude for m = +1 and 
m - +{. Third row: Amplitude difference between retrograde 
and prograde waves of the same azimuthal number. Fourth row: 
Net luminosity for a given azimuthal number. 



this shear layer is crucial in filtering the low order waves that 
travel through it. The net wave flux across the shear layer de- 
pends of the asymmetry of the filter and thus, on the difference 
between the rotation velocity below the shear layer and that of 
the convection zone. Let us note however that the thickness of 
the oscillating layer is independent of this asymmetry. Last but 
not least, wave momentum transport does not have a diffusive 
behavior and as such must not be treated merely as a turbulent 
diffusion process. 

Figure |6] illustrates differential damping in a shear region, 
with both a "rapid" and a "slow" layer. The m — mode is not 
affected by the differential rotation 5 , and, for a given frequency, 
the damping length 6 is inversely proportional to 1(1 + 1). The 
"rapid" layer increases the damping of prograde modes, and 
reduces it for the retrograde modes. The reverse is obtained 
from the "slow" layer. 

FigureQshows the local momentum deposition for various 
values of differential rotation. As far as differential rotation is 
not too large, momentum deposition varies linearly with <5D. 
Actual calculations of the evolution of the distribution of an- 
gular momentum show that, for realistic values of braking (ac- 
cording to Kawaler 1988), the star remains in this linear regime. 



This mode does not carry angular momentum. 
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Fig. 6. (top) Instantaneous rotation profile (middle, bottom) 
Instantaneous damping length (corresponding to a reduction of 
wave amplitude A = Ao/e) for wave frequencies of v = 0.5 and 
v = 1 yuHz in a 1.3 M G star. The m = mode (open dots) is 
not affected by rotation, the m-+( (squares and triangles) are 
the most affected. Below the shear layer, the radiative region is 
rotating more rapidly; it is the retrograde modes of low order 
that penetrate deeper, leading to a net deposition of negative 
momentum in the interior. The low (I — 1 and 2) modes have 
to travel back and forth a few times before they are damped, 
and are not shown in the bottom graph. 

3.3. Energy Considerations 

It has been shown elsewhere (Kumar & Quataert 1997, Zahn 
et al. 1997, TKZ) that gravity waves can carry enough angular 
momentum to slow down the radiative zone of low mass stars 
on time-scales of order 10 7 years. One may wonder however if 
the deposition of energy by waves may have an impact on the 
stellar structure. Let us first begin by comparing the amount of 
energy contained in gravity waves with other quantities. This 
will be done here for a ZAMS 1 M model. 

The total energy luminosity in waves is of order 



10 29 



ergs s , while the total luminosity at the base of the con- 



vection zone is 2.5 x 10 33 ergs s , of which a large part is con- 
vective (see Fig . |3} - The wave luminosity thus represents about 
0.01% of the convective luminosity. 

Let us next look at the amount of kinetic energy that is 
stored in rotation Kq. We have Kq = jIQ? with / the moment 
of inertia of the star. For the whole star, it is of order 10 54 cm 2 g. 
For an initial rotation velocity of 100 km s _1 , this leads to an en- 
ergy of 10 46 ergs. As the star is spun down, this energy must be 
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Fig. 7. Local deposition of angular momentum for various val- 
ues of differential rotation. The rotation profile rises linearly 
between the base of the convection zone and (l-m/M) = -2.85 
by a value of 5CI, which varies logarithmically from 10~ 5 to 
10 fiHz. This is illustrated for a 1.3 M Q model. For moderate 
amounts of differential rotation (below ~ 10~ 2 //Hz), the net 
momentum flux at a given point is a linear function of differen- 
tial rotation. 



dissipated and will be added to the thermal energy. However, 
the total rotation energy is two order of magnitude lower than 
the thermal content and so, the impact will be negligible on the 
stellar structure. 

4. A Model for Angular Momentum Evolution by 
Waves 

Gravity waves lead to two different features that must be incor- 
porated in stellar evolution codes: 

- They produce a shear layer, that generates turbulence close 
to the bottom of the convection zone; 

- They deposit negative (positive) momentum throughout the 
radiative interior when the convection zone rotates slower 
(faster) than the radiative zone. 

The turbulent region is always present, and the shear layer 
oscillation is (almost) independent of the presence or absence 
of differential rotation below. The magnitude of the turbulence 
and the size of the turbulent region are self-regulated and de- 
pend on the wave flux. For a given wave-excitation model, and 
a given model for the turbulent diffusion, we obtain a localized 
turbulent region. 

The net momentum luminosity spectrum below the shear 
layer for various masses is shown in Fig. [8] To obtain the net 
momentum deposition, one must follow the local momentum 



luminosity 



-Cj(r) = ^ £j(, m (fshear layer) exp [-r(r, 0", £)] 



(6) 



cr,Lm 



8 



S. Talon & C. Charbonnel: Stellar models including rotation and IGW I. 



and where the local damping rate takes into account the mean 
molecular weight stratification 



T(r,crJ) = [€(€+1)] 



where A^ 2 



(K + v,) 



AW 2 



a 4 \N 2 -cr 2 



dr 



(7) 



Nj + N 2 is the total Brunt-Vaisala frequency, Nj 
is its thermal part and N 2 is due to the mean molecular weight 
stratification (cf. Zahn et al. 1997). <x is the local, Doppler 
shifted frequency 



<r(r) — a) - m (£l(r) - O cz ) 



(8) 



and a) is the wave frequency in the reference frame of the con- 
vection zone. Only those waves with a sufficient amplitude af- 
ter filtering have to be traced. 

When meridional circulation, turbulence and waves are 
taken into account, the evolution of angular momentum thus 
follows 



[pr A Q.u\ 



J_d_ 

5r 2 dr 



3 1 8 r m 
8n r z dr 



]_d_ 

r 2 dr 



p(Vt + Vwaves)?" 



,5Q 



(9) 



here p is the density, U the radial meridional circulation veloc- 
ity, v, the turbulent viscosity due to differential rotation away 
from the shear layer, and v waV es the diffusion coefficient as- 
sociated with wave-induced turbulence (see Eq. [n]in § 7.1). 
Horizontal averaging has been performed for this equation, 
and meridional circulation is considered only at first order (see 

§ED- 



5. Evolution of Chemicals 

In most hydrodynamical models, chemicals evolve under the 
action of meridional circulation, turbulence and atomic diffu- 
sion. All these topics are discussed elsewhere (see e.g. Talon 
2004 and references therein for a discussion of rotational mix- 
ing, and Turcotte et al. 1998 and references therein for a review 
of microscopic diffusion processes). However, when including 
gravity waves in momentum transport, one must also include 
their direct contribution to the transport of chemical species 6 . 
This subject has received attention from several authors and we 
will here describe the main mechanisms involved. 

The first process we shall discuss here is the turbulence 
generated in the shear layer by the shear instability. In a first 
step, energy is transfered by waves from the convection zone 
to the shear layer and stored into differential rotation, which 
can then be converted to turbulence by the shear instability. 
This view is similar to the idea of Canuto (2002) that gravity 
waves act as a source term in the equation that describes tur- 
bulence. This is also similar to the mixing described by Young 
et al. (2003). In our framework, turbulent diffusion is evaluated 
by taking the average of Eq. over several SLO. 

Weak mixing can also be induced by second order mass- 
transport effects in a diffusive medium (Press 1981, Schatzman 



6 This adds to the indirect effect described above, that is the modi- 
fication of the rotation profile which changes rotational mixing. 



1993, Montalban 1994). A diffusion coefficient can then be as- 
sociated with wave dissipation by combining the average wave 
velocity and the average damping length of waves. If the con- 
comitant transport of angular momentum is taken into account, 
the damping length of waves is reduced and the effect is to 
somewhat lower the size of the region over which this process 
is efficient. The thickness of the shear layer should, in many 
instances replace the damping length calculated in the case of 
solid body rotation. Differential rotation thus reduces the extent 
of this effect; its magnitude is at maximum of the same order 
as turbulence induced by the SLO. Considering the uncertain- 
ties in waves fluxes as well as on the mixing induced by those 
motions we suggest to neglect this effect altogether. 

Finally, in certain circumstances, wave amplitude can rise 
to the point of becoming non-linear and inducing shear-mixing 
(Press 1981, Garcia Lopez & Spruit 1991). The non-linearity is 
given by 



e = k,£ h = 



(10) 



(Press 1981). The wave becomes non-linear when this param- 
eter is of order 1 or more. In terms of energy luminosity, it is 
given by 



4npr 2 (JV 2 - o- 2 ) 



2\i/2 



1/2 



(ID 



As estimated by Press (1981), some waves are slightly non- 
linear just below the convection zone; this non-linearity re- 
mains only over one damping length. However, since this re- 
gion is already highly turbulent, this small non-linearity does 
not greatly modifies wave damping. Close to wave turning 
points (where cr — > 0) however, non-linearity may again be- 
come important. This will affect mostly low degree, low fre- 
quency waves by increasing their damping in that region. 
However, since radiative damping is already strong there, it will 
not greatly modify wave amplitudes. However, in some cases 
it could increase the amount of mixing of chemicals. Garcia 
Lopez & Spruit (1991) showed this has a significant effect only 
in a small part of the HR diagram (namely in stars close to the 
Li dip for a specific choice of the mixing length parameter). We 
will discuss the relevance of this effect in a future paper, where 
models of various masses will be implemented. 



6. Waves in a Static Model 

Before we apply our filtered wave model to an evolutionary 
calculation, let us begin by discussing some more the results 
obtained by TKZ in the case of a static model. Those are based 
on a rather small differential rotation, and use a given wave 
spectrum. It is worth verifying how sensitive are those results 
to other initial profiles or different distributions of energy. In the 
models presented in this section, we show results obtained for 
small time steps (of 1 year, as in TKZ), following the details of 
the shear layer oscillation. The only other process for momen- 
tum transport that is considered here is shear turbulence, and 
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Fig. 8. Filtered momentum luminosity spectrum for the m = ( mode (that is Xnet = -C m =-c - -C m =c) below the shear layer for 
various masses. Differential rotation is 6Q. = 1CT 4 /iHz over 5% in radius. 




Fig. 9. Evolution of the rotation profile in the case of very 
strong and moderate initial differential rotation in a 1 M star. 

no surface braking is applied 7 . Contrary to the TKZ study, here 
turbulence is directly related to the local shear rate (cf. Eq.|5jl. 



x ~ 
id 



initial profile 
1 000 yrs 
20 000 yrs 




Fig. 10. Evolution of the rotation profile with moderate initial 
differential rotation for a wave spectrum derived from GMK 
spectrum (cf. Eq.[0 and for a flat spectrum with the same total 
energy flux. Calculations are made in a 1 M star 



6.1. Strong Initial Differential Rotation 

Let us discuss some more results presented in § I3.2I In Fig. [9] 
we compare the evolution of the rotation profile in cases of very 
strong and moderate 8 differential rotation. Momentum extrac- 
tion from the core is clearly visible in both case. The core's 
slow rotation slowly propagates toward the surface. It does so 
with a velocity that gets smaller as the front progresses to a 
region where the local angular momentum (oc r 2 ) is larger. 

7 However, the initial rotation profile has a surface convection zone 
rotating slower than the radiative interior. 

8 This is the same as was used by TKZ. 



Differential rotation at the core boundary remains larger than 
in TKZ because the local turbulence v, is smaller 9 . The impor- 
tant point to note is that in the presence of a large differential 
rotation, the local frequency of retrograde waves becomes very 
large and their damping is largely reduced in the inner regions. 
However, once the core has been spun down which is easily 
done since it contains very little momentum, a "slow" front can 
propagate toward the surface. 



9 In TKZ, the local turbulence was fixed (as in Eq.[4j and did not 
depend on the local shear rate. 
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6.2. Wave Spectrum 

Another delicate point is the form of the wave spectrum. As 
already mentioned (§|2j, the exact spectrum produced by tur- 
bulence in the convection zone is not well constrained, and 
various analytical and numerical studies lead to different pre- 
scriptions. It is thus important to understand how modifying 
the spectrum influences the results. Figure ITU1 illustrates this 
comparison. The flat spectrum produces a front that migrates 
somewhat faster, but deposits less momentum in the interme- 
diate region. In that case, the shear layer does not oscillate, 
but still produces a filter that preferentially damps prograde 
waves. While the exact spectrum changes the results quanti- 
tatively, qualitatively global effects are similar. 

7. Waves in an Evolutionary 1.2 M Q Model 

7.1. Code Description 

General inputs for stellar physics 

Stellar models are computed with the stellar evolution code 
STAREVOL (Forestini 1991; Siess et al. 1997, 2000; Palacios 
et al. 2003). Our equation of state follows the Pols et al. 
(1995) formalism. Thermodynamical features of each plasma 
component (ions, electrons, photons and H2) are obtained by 
minimizing the Helmholtz free energy that includes separately 
non-ideal effects, and allows to treat ionization analytically, 
leading to smooth profiles for thermodynamical quantities. 
Radiative opacities are taken from Alexander & Fergusson 
(1994) below 8000 K and from Iglesias & Rogers (1996) at 
higher temperatures. Our nuclear reaction network follows 
53 species (from 'H to 37 C1) through 180 reactions. Nuclear 
reaction rates have been updated with the NACRE compilation 
(Angulo et al. 1999). Convection is treated according to 
the mixing length formalism with a p = l/H p = 1.75. No 
overshooting is considered. 

Rotation 

We follow the evolution of the rotation profile from the zero 
age main sequence on, assuming initial solid body rotation. 
The surface rotation velocity on the ZAMS is taken equal to 
50 kms _I . Surface spin-down follows Kawaler (1988) includ- 
ing saturation at f2 sat = 1.5x 10~ 5 rad s _1 (or v ~ 10 km s _1 ) 



dt 



-K s Q?(RIR & ) lll (MIM Q Y 112 



(Q < Q sat ) 



K s nn z s , t (R/R ) l/2 (M/M Q r llz (Q > Q sat ) 



-1/2 



(12) 



The issue of saturation in the context of meridional circulation 
is discussed at length in Palacios et al. (2003). 

IGW 

To incorporate IGW, we first calculate filtered lu- 
minosities (cf. Fig. |HJ for fixed differential rotations 
(6D. = 10~ 4 , 10~ 3 , 10~ 2 and 10~' fMz) in static ZAMS 
models 10 . The wave luminosity below the surface convection 

10 We showed in Talon & Charbonnel (2003) that in the type of star 
we consider for this full evolutionary computation (i.e., Pop I main se- 
quence star with a mass of 1.2M Q ) the depth and structure of the con- 



zone is then linearly interpolated from those tables for the 
actual differential rotation just below the shear layer. The local 
momentum luminosity is obtained by calculating the damping 
integral for each individual wave (cf. Eq.[7Ji and then summing 
over all waves (cf. Eq.[6}. 

Meridional circulation 

Meridional circulation is treated as an advection process for the 
transport of angular momentum, assuming strong horizontal 
turbulence and follows Zahn (1992). 

Turbulence and diffusion 

Turbulence is assumed to be strongly anisotropic, and every- 
where we will assume that turbulent diffusion is equal to the 
turbulent viscosity. The vertical component of the turbulent vis- 
cosity 

« fl/ crit (rdQ/dr) 2 



D\ = v v = - 



5N 2 T /(K + D h )+N 2 /D h 



(13) 



takes into account the weakening effect of thermal diffusivity 
(K) on the thermal stratification and of horizontal turbulence 
(Z?/,) on both the thermal and mean molecular weight stratifica- 
tions (Talon & Zahn 1997). 
Horizontal turbulence follows Zahn (1992) 

d(pr 2 U) U din r 2 Q. 



D, 



Vh 



r 



1 



3pr dr 



2 dlnr 



(14) 



with C/, = 1. The evolution of momentum is then given by 
Eq. ©. 

Transport of chemicals 

For chemicals, the combination of meridional circulation and 
horizontal turbulence results in a vertical effective diffusion 



£>eff = 



\rU(r)\ z 
30 D h 



(15) 



(Chaboyer & Zahn, 1992). 

Atomic diffusion driven by gravitational settling and ther- 
mal gradients is included using the formulation of Paquette 
etal. (1986). 

The evolution of an element c, then follows 



Ac, 



1 d 



dt 



dr 



r 2 p I^UipCi + (D eff + D v + D waves ) ^ 



(16) 



with c,- the nuclear production/destruction rate and f/, ; , the mi- 
croscopic diffusion velocity with respect to protons. 
Wave induced turbulence is treated explicitly through coeffi- 
cient D waves that we assume equal to the coefficient v waves which 
enters Eq.|5] To evaluate this coefficient, we perform a time av- 
erage of the diffusion coefficient associated with the local shear 
instability (Eq. I13l > over several SLOs (see Fig.@J and § I3.1> . 
For evolutionary calculations, we use an analytical fit to this 
curve. In this application to a 1.2 M star, we used 

log D waves = (log K + 2 lo g L - 68) exp [(r cz - rf /o] (17) 

vective envelope, and thus the wave characteristics, do not vary sig- 
nificantly over a main sequence lifetime. It is thus justified to perform 
the first test with the spectrum characteristics of the ZAMS models 
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Fig. 11. Evolution of the rotation profile in a complete model 
where the transport of angular momentum is due to inter- 
nal gravity waves, meridional circulation and turbulence. The 
model shown is for a 1 .2 M Q Z = 0.02, star with an initial rota- 
tion velocity of 50 kms -1 . The curves are labeled according to 
the corresponding ages in Gyr. 

with 

cr = 0.0003 (1 - r cz ) . (18) 

This fit is compared to the time-average in Fig.|4] Let us remark 
that above ~ 10 8 cm 2 s~', mixing may be considered "instan- 
taneous" compared with evolutionary time-scales and thus, our 
analytical fit is aimed at reproducing the numerical turbulent 
viscosity below this value. 

7.2. Evolution of the Rotation Profile 

Let us first concentrate on the evolution of the rotation profile 
when momentum deposition by IGW is taken into account in 
conjunction with shear turbulence and meridional circulation. 
Here, the short time-scale SLO is present only as a filter; the 
magnitude of the wave flux depends on the amount of differ- 
ential rotation between the base of the convection zone and a 
region just below the SLO (see § I3.2> . In quasi-solid body rota- 
tion (with the surface rotating slightly slower because of brak- 
ing), low degree waves penetrate all the way to the stellar core, 
and are damped (and thus deposit their negative momentum) 
over the whole radiative region. However, since the amount of 
angular momentum contained in the radiative core is minute, 
the local deposition of even a small amount of momentum is 
enough to spin it down significantly. In a "slow" region, damp- 
ing of retrograde waves increases (cf. Eqs. 0and[S}, and this 
leads to the formation of a "front", which propagates from the 
core to the surface. The propagation of a first front is seen in 
Fig. CD (curves at 0.2, 0.25, 0.3 and 0.35 Gyr). A second front 
evolves from 0.4 to the 0.7 Gyr curves. 
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Fig. 12. Same as Fig.[^f or 0.2, 0.5 and 0.7 Gyr in a case with- 
out (left) and with (right) internal gravity waves included. The 
time sequence goes from the highest to the smallest surface ve- 
locity 

At the age of the Hyades, differential rotation is consider- 
ably reduced. In particular, it is of interest to compare calcu- 
lations including and not including waves. In the second case 
(left plot in Fig. I12> . the amount of differential rotation at the 
age of the Hyades is very large. This is in in agreement with the 
results obtained in the solar case by Matias & Zahn (1998) who 
performed calculations under the same hypothesis as here, and 
by Chaboyer et al. (1995) who approximated meridional circu- 
lation as a diffusive process. 

The present complete model confirms the ability of gravity 
waves to efficiently extract angular momentum from the deep 
interior of solar-type stars. It shows how the momentum redis- 
tribution proceeds when the stars are spun down via magnetic 
torquing. 

Let us see now what are the consequences on meridional 
circulation and shear turbulence, and on the Li depletion due to 
rotational mixing. 

7.3. Meridional Circulation Velocity and Diffusion 
Coefficients 

Figure [D] presents profiles of the vertical component of the 
meridional velocity U at 0.5 Gyr in the 1.2 M models with 
and without gravity waves. When not including waves, there 
are two circulation loops. The meridional velocity is negative 
in the external part of the radiative zone down to ~ 0.3 M„. This 
meridian loop brings matter upward at the equator and down 
in the polar regions, in response to the extraction of angular 
momentum due to braking. Deeper the circulation is positive 
(bold line), indicating an inward transport of angular momen- 
tum. When internal waves are taken into account, several loops 
of circulation appear, with negative and positive loops (light 
and bold lines respectively) alternating. 

These differences reflect directly in the transport coeffi- 
cients which are shown at the same age (0.5 Gyr) in Figs. [21 
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Fig. 13. Profiles of the vertical component of the meridional 
velocity in a 1.2 M Q star at 0.5 Gyr in a simple rotating model 
(dotted line) and in a rotating model including gravity waves 
(full line). The bold lines indicate positive values of U . 



and^] Details in the profiles vary with time but this is a typ- 
ical illustration. Strong variations (i.e. peaks and gaps) in the 
profiles of both effective and total diffusion coefficients reflect 
the presence of several circulation loops, each drop correspond- 
ing to an inversion of the flow direction. When gravity waves 
are included, the amplitude of these coefficients is reduced and 
vertical turbulence (D v ) is less developed than in the pure ro- 
tating models because of the decrease of the overall differential 
rotation. 

In Fig. the wave-induced turbulence (Z) waves ) used in 
the present simulation with IGW and given by Eq. El is a l so 
shown. As can be seen, this coefficient drops very rapidly be- 
low the convection envelope and is much smaller than the to- 
tal diffusion coefficient coming from rotation except very close 
from the convection region. 



Rotation 
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Fig. 14. Comparison of the transport coefficients in a 1 .2 M 
star at 0.5 Gyr in a simple rotating model (left) and in a rotating 
model including IGW (right). Coefficients are identified on the 
figure. D v is the dashed line and D e jf a full line. 




7.4. Evolution of Chemicals 

Our goal in the present paper is not to make detailed compari- 
son of the model predictions with all the observable constraints 
that come from abundance anomalies in stars. This will be done 
in forthcoming studies. We wish however to illustrate briefly 
the behavior of 4 He and 7 Li in our 1 .2 M Q rotating model in- 
cluding gravity waves. The case of 4 He is related to possible 
impact on the overall evolution and lifetimes, and in addition 
it illustrates the interaction with atomic diffusion. On the other 
hand 7 Li is a fragile element which helps to probe the status of 
the external stellar layers. 

The 4 He profile at the age of the Hyades is shown in Fig. 1161 
for both the complete model and the model without IGW. One 
sees there two effects. First, at this age, the model with IGW 
is slightly less evolved (i.e. it has a lower central 4 He content). 



Fig. 15. Total diffusion coefficient in a 1 .2 M e star at 0.5 Gyr 
in a simple rotating model (dotted line) and in a rotating model 
including IGW (full line). The dashed-dotted line gives £> W aves- 



Second, its He surface abundance is slightly lower and the He 
gradient just below the convection envelope is slightly steeper. 
This is due to the lower amplitude of the total diffusion coeffi- 
cient in the models including waves, which allows 4 He to settle 
more under the effect of atomic diffusion. One sees here that 
Owaves has a negligible effect on the 4 He behavior and does not 
compensate for the decrease in the total diffusion coefficient in 
the external part of the radiative zone. 

The best observational constraints available to test our pre- 
dictions are the Li data in open clusters. As discussed in § 1, 
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Fig. 16. 4 He profile at the age of the Hyades in the complete 
model (full line) and in the model without the IGW (dashed 
line) for a 1 .2 M Q star. 

models in which the transport of angular momentum is carried 
out only by meridional circulation and turbulence fail to repro- 
duce the rise of the Li abundance on the right side of the Li dip. 
This is confirmed in the present model without IGW which lies 
in this region and has an effective temperature of r e ff = 6220 K 
at the age of the Hyades. As can be seen in Fig. El the 7 Li 
surface abundance in this model (dotted line) is ~ one order of 
magnitude smaller than the Li value in stars of 6220 ± 100 K 
in the Hyades and which is indicated by the box. For the model 
including IGW (full line), the magnitude of both meridional 
circulation and turbulence is reduced. Consequently 7 Li is less 
destroyed and our prediction accounts nicely for the data. Note 
that the surface lithium decrease in these stars is not dominated 
by atomic diffusion, but still by the (reduced, compared to the 
case without IGW) rotational transport of this element down 
to regions where it is nuclearly destroyed. Again, the effect of 
£>waves is negligible. 

In a forthcoming paper, we will present models of Pop I 
stars with various initial masses and initial rotation rates and 
compare in details the Li predictions on the red side of the Li 
dip with observations in open cluster and field stars. 

8. Discussion and Conclusions 

In this paper, we examined how it is possible to include all the 
effects of internal gravity waves into evolutionary calculations. 
The main challenge is that IGW tend to produce a thin shear 
layer that oscillates with a very small time-scale. Talon et al. 
(2002) showed that this shear layer oscillation (or SLO) acts 
as a filter for waves, being more transparent to waves that will 
reduce the differential rotation between the convection and the 
radiation zones. We verified that this filter is linear in SO. and as 
such, details of the SLO need not be considered for variations 
occurring over long time-scales. 
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Fig. 17. 7 Li profile at the age of the Hyades in the complete 
model (full line) and in the model without the IGW (dashed 
line) for a 1.2 M star. The box indicates the observed Li 
abundances in Hyades stars with an effective temperature (± 
100 K) corresponding to that of the model at 0.7 Gyr (data by 
Boesgaard & King 2002) 



We also presented the first evolution model that includes 
the hydrodynamical processes induced by rotation and internal 
gravity waves. We focused on a 1.2 M , Z = 0.02 star which 
lies on the red side of the Li dip (T e g = 6220 K at 0.7 Gyr). For 
the rotation profile, IGW lead to the appearance of a "slow" 
front, propagating from the core to the surface. When the outer 
convection zone is constantly spun down, several fronts prop- 
agate; this propagation is very rapid at the beginning because 
differential rotation is large, and slows down with time. 

For the evolution of chemicals, Talon & Charbonnel (1998) 
expected IGW to reduce rotational mixing together with differ- 
ential rotation. This effect is confirmed in this first fully con- 
sistent study. The surface of the Li abundance in our complete 
model is in perfect agreement with the data in the Hyades stars 
of similar effective temperature. 

Let us remind to the reader that, while the formalism pre- 
sented in this paper is quite general and allows to properly in- 
clude IGW in complete hydrodynamical stellar models, sev- 
eral uncertainties remain. Firstly, a very delicate issue is that 
of IGW generation. The GMK model most certainly underes- 
timates the wave flux since it considers only bulk excitation, 
and not overshooting or convective penetration. However, since 
wave generation has to be proportional to the convective lu- 
minosity, differential properties between different stellar types 
should vary less. We thus expect the wave flux to be larger 
than calculated, but by a similar amount in each stellar type. 
Complete studies for different stellar types should permit to 
calibrate this amount, that could then be used uniformly. 

The second uncertainty is the exact value for the turbulent 
viscosity to use in the SLO (v waves , Eq.EJ- Here, we estimated 
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a viscosity based on the shear instability and assumed a small 
amount of enlargement and averaging in order for the turbulent 
viscosity not to be too local. However, this model has only a 
rather small impact on numerical results since in all cases it 
remains very localized below the convection zone. 

The last point to bear in mind is that no latitudinal depen- 
dence is considered here; only a horizontal average is used. Let 
us remark that, in order to be able to take into account this de- 
pendence, observational data of latitudinal differential rotation 
in stars other than the Sun must be obtained. 

The results on the evolution of the rotation profile and on 
the lithium abundance presented in this paper are very encour- 
aging and should now to be confronted to the data for stars 
of various masses, metallicities and ages. This will be done in 
forthcoming studies. 
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